{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pulsar Gating Tutorial\n",
    "\n",
    "You may be looking at a supernova remnant or seeking new pulsar wind nebulae and you might need to remove the accompanying pulsar. Or, you could be looking at starburst galaxies at higher galactic latitudes, and a pesky millisecond pulsar keeps getting in your way.\n",
    "\n",
    "Pulsar gating is the process of removing the flux of a bright pulsar from a region of interest, thus enabling the analysis of fainter sources in the region. You can accomplish this by assigning pulse phases to the pulsar by using a known ephemeris. You will also need to use the pulsar light curve to define an off-pulse interval in phase-space. Once the phases are assigned, you can use **gtselect** to filter out the unwanted phase periods, effectively removing the pulsar from the data completely."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Step 1: Extract your data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this tutorial, we assume you know how to extract your data, prepare your data for analysis, and explore your data.\n",
    "\n",
    "First, let's take a look at a set of data that contains a bright pulsar.\n",
    "\n",
    "We'll use the pulsar Geminga (98.476204, 17.770661) with a 15° ROI, and the time range from the START of the mission to MJD 55196 (these position and time range are designed to match the data in the ephemeris file used in the next step).\n",
    "\n",
    "So, for the search, you will put in:\n",
    "- Geminga (98.476204, 17.770661) (J200)\n",
    "- 15° ROI\n",
    "- Observation dates: (Start, 55196) MJD\n",
    "- Energy Range (100, 500000)\n",
    "\n",
    "Geminga is very bright in the LAT, so you will get several event files to cover the full time span.\n",
    "Below, we extract the Geminga data from the [LAT Data Server](http://fermi.gsfc.nasa.gov/cgi-bin/ssc/LAT/LATDataQuery.cgi)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src='https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/images/Pulsars/Geminga_query.png'>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Download the files into your current directory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!mkdir data\n",
    "!mv *.fits ./data\n",
    "!mv ./data/*_SC00.fits ./data/spacecraft.fits"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will also want a list of event files:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!ls ./data/*PH*.fits > ./data/Geminga.txt\n",
    "!cat ./data/Geminga.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we want to use [gtselect](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtselect.txt) to merge the event files and make the zenith cut.\n",
    "\n",
    "The data server already trims 30 seconds off each end of the LAT data to ensure that the spacecraft file will contain entries beyond both ends of the data, which is necessary for proper use of the pulsar tools. You can find the start and end times of the full dataset by looking at the TSTART and TSTOP keywords in the header of the spacecraft file."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src='https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/images/Pulsars/Geminga_header.png'>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run **gtselect**. Here, `Geminga.txt` is a listing of all the event files to allow them to be merged. The spacecraft file (*SC*) has been renamed to `spacecraft.fits` for simplicity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "gtselect evclass=128 evtype=3\n",
    "    @./data/Geminga.txt\n",
    "    ./data/Geminga.fits\n",
    "    INDEF\n",
    "    INDEF\n",
    "    INDEF\n",
    "    INDEF\n",
    "    INDEF\n",
    "    100\n",
    "    500000\n",
    "    90"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then run **gtmktime** to update the GTIs, otherwise the exposure calculation will be off.\n",
    "\n",
    "While **gtmktime** can compensate for some exposure changes, it does not handle phase cuts. So you will need to scale any outputs from the analysis by hand. We'll talk about how to do that at the end of this tutorial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "gtmktime\n",
    "    ./data/spacecraft.fits\n",
    "    (DATA_QUAL>0)&&(LAT_CONFIG==1)\n",
    "    no\n",
    "    ./data/Geminga.fits\n",
    "    ./data/Geminga_gtis.fits"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, use **gtbin** followed by *ds9* to take a look at the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "gtbin\n",
    "    CMAP\n",
    "    ./data/Geminga_gtis.fits\n",
    "    ./data/Geminga_cmap.fits\n",
    "    ./data/spacecraft.fits\n",
    "    300\n",
    "    300\n",
    "    0.1\n",
    "    CEL\n",
    "    98.4756\n",
    "    17.7703\n",
    "    0\n",
    "    CAR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!ds9 ./data/Geminga_cmap.fits"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is the counts map showing a very bright Geminga centered in the ROI.\n",
    "\n",
    "<img src='https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/images/Pulsars/Geminga_On.png'>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Step 2: Finding the pulsar ephemeris\n",
    "\n",
    "Now you have data with a very bright pulsar in it, and you need to remove this pulsar from your data.\n",
    "\n",
    "To do so, you have to find out a few things about that pulsar. This information can be found in the pulsar's ephemeris file and in any accompanying paper, both of which can be found at the [LAT Data Access → LAT Pulsar Ephemeris](http://fermi.gsfc.nasa.gov/ssc/data/access/lat/ephems/) page:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import Image,HTML\n",
    "display(HTML(\"<iframe src='http://fermi.gsfc.nasa.gov/ssc/data/access/lat/ephems/' width='1000' height='500'></iframe>\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This page contains the ephemerides that have been used by the LAT team for their published papers to date. It provides several useful bits of information: the ephemeris file for each pulsar, the timing contact (i.e. the person who actually calculated the ephemeris), links to each paper (papers are likely to contain light curves for each pulsar), and finally the contact authors for each paper.\n",
    "\n",
    "If you have questions about the ephemeris itself, you should contact the timing contact. For question about light curves and interpretation of the pulsar data, contact the paper's contact authors. All email addresses use \"at\" for the @ symbol, and \"dot\" for the periods, just to keep email addresses less obvious on the open internet. You'll have to correct those before you send any emails.\n",
    "\n",
    "The ephemeris files themselves are listed in the \"Table of Published Ephemerides\" under the \"Timing Files\" heading. There are two types of ephemeris files: `.par` files, which are used for the TEMPO2 tool, and `D4` files which are used by the Fermitools. As a note, the TEMPO2 tool is not associated with the FSSC. However, members of the LAT instrument team have provided specialized files to allow LAT data to interface with that toolset."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's download an ephemeris file and see what it looks like. Scroll to the first row of the table that contains a D4 in the second column. This should be the pulsar Geminga, a bright gamma-ray pulsar that is radio quiet."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src='https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/images/Pulsars/Geminga_ephem.png'>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the purposes of this tutorial, these files can also be downloaded below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!wget https://fermi.gsfc.nasa.gov/ssc/data/access/lat/ephems/0633+1746/0633+1746_ApJ_720_272_2010_D4.fits\n",
    "!wget https://fermi.gsfc.nasa.gov/ssc/data/access/lat/ephems/0633+1746/0633+1746_ApJ_720_272_2010.par"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!mv *.fits *.par ./data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you open the `.par` file associated with Geminga using a text viewer, you will see some important information:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. **Pulsar name and position** — these positions are typically determined by the timing data, but sometimes come from X-ray or radio interferometric observations. For some pulsars there may be proper motion (PMRA, PMDEC) or parallax (PX) values as well.\n",
    "\n",
    "\n",
    "2. **Spin information** — this gives the primary spin frequency (F0), and one or more derivatives (note that the spin frequency is equal to F0 at PEPOCH). The first derivative (F1) is the rate at which the pulsar is slowing down and is dominated by the secular spindown of the pulsar. Any higher order derivatives are usually dominated by timing noise. If the pulsar you are looking at is in a binary, you will also see orbital parameters in this section. Some par files include \"WAVE\" terms, which are used to model timing noise using harmonically related sinusoids (Hobbs et al. 2004).\n",
    "\n",
    "\n",
    "3. **START and FINISH** — These provide information about the time period over which the ephemeris is good. Some pulsars (particularly millisecond pulsars) are very stable, and an ephemeris can be safely extrapolated for months beyond the FINISH time. However, many gamma-ray pulsars have significant timing noise or glitches, which make it very unsafe to extrapolate. If you see frequency derivatives of F2 or higher or WAVE parameters in the PAR file, DO NOT use it beyond its valid dates! If you cannot find an ephemeris that covers the time span you need, email the timing contact to find out if there is something more recent.\n",
    "\n",
    "\n",
    "4. The `TZRMJD`, `TZRFRQ`, and `TZRSITE` parameters define phase 0.0 for the ephemeris. The fiducial point on the radio profile arrived at `TZRSITE` at frequency `TZRFRQ` at the moment `TZRMJD`. Often, but NOT always, the fiducial point on the radio profile is the peak, but you must always ask the person who made the ephemeris to be sure, because other conventions are used!\n",
    "    * WARNING: If your PAR file comes from Tempo2 (usually you can tell by \"EPHVER 5\" in the file), then the default time system is TCB, NOT TDB! Tempo2 only uses TDB if there is \"UNITS TDB\" in the par file. This is VERY important if you feed the values in the par file to a program that expects TDB, such as tempo1 or the D4 FITS files."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Step 3: Assigning pulse phases\n",
    "\n",
    "Now that you have the ephemeris, you need to use this information to assign pulse phases to the events in your data file.\n",
    "\n",
    "Run [gtpphase](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtpphase.txt) to add a column containing pulse phases to your events file. It may be advisable to save a backup copy of your events file before running this tool as it permanently modifies the input file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!cp ./data/Geminga_gtis.fits ./data/Geminga_gtis.fits.backup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "gtpphase\n",
    "    ./data/Geminga_gtis.fits\n",
    "    ./data/spacecraft.fits\n",
    "    ./data/0633+1746_ApJ_720_272_2010_D4.fits\n",
    "    J0633+1746\n",
    "    DB"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The pulsar name needs to match the value in first column of the D4 ephemeris file. Here, Geminga is actually called \"J0633+1746\" in the file.\n",
    "\n",
    "<img src='https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/images/Pulsars/Geminga_D4.png'>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The file `Geminga_gtis.fits` will now have a new column, PULSE_PHASE, which you can use to filter your data. But it's always wise to check, so verify the data looks good by creating a quick folded light curve using the FV histogram function.\n",
    "\n",
    "<img src='https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/images/Pulsars/Geminga_phaseogram.png'>\n",
    "\n",
    "Here I have used a bin size of 0.005 in phase space. You should see two strong peaks from the pulsar, with noise between the peaks. That noise is the data we're interested in. The peaks are the pulsed flux that we wish to remove, which we will discuss in the next step. If you don't see this pattern, you may have made a mistake in assigning phases."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you want to assign pulsar phases with the Fermitools but don't have a D4 database file, you can use [gtpphase](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtpphase.txt) and enter the timing parameters from the .par file by using the \"FREQ\" option."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!cp ./data/Geminga_gtis.fits.backup ./data/Geminga_gtis_par.fits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "gtpphase\n",
    "    ./data/Geminga_gtis_par.fits\n",
    "    ./data/spacecraft.fits\n",
    "    none\n",
    "    J0633+1746\n",
    "    FREQ\n",
    "    54800\n",
    "    MJD\n",
    "    TDB\n",
    "    98.476204\n",
    "    17.770661\n",
    "    0\n",
    "    4.2175670649262114748\n",
    "    -1.9525033316064710553e-13\n",
    "    0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, the Fermitools do not support phase assignment for pulsars with more than one WAVE term. So if your .par file has timing parameters at higher orders than F2, you will need to use TEMPO2 to assign pulse phases. A discussion of TEMPO2 phase assignment is available in the second section of this tutorial."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Step 4: Remove events from the pulsar"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you look closely at the light curve, you will see that during the period when the pulsar is not pulsing, there are still events coming in.\n",
    "\n",
    "These \"off-pulse\" events are the data you care about, so we want to filter out the pulsar \"on-pulse\" events.\n",
    "\n",
    "To do this, you define an \"off-pulse\" interval. The off-pulse could be both phase intervals between the pulses (in this case, .22-.58 and .72-.05). But some pulsars never turn off completely between the pulses. In fact, for Geminga you'll notice the background level after the second pulse is quite a bit higher than between the pulses. This higher level is called \"bridge emission\" and must also be removed in order to completely remove the pulsar.\n",
    "\n",
    "For pulsars with an interpulse, you will probably only be able to define a single off-pulse interval."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You might want to look at the LAT team paper at this point. For many of the sources, the people who did the pulsar analysis also defined on off-pulse period for the pulsar. In either case, pick what you think is the best off-pulse interval and move on to the next step. For this analysis we will use phases .22-.58 for our off-pulse period.\n",
    "\n",
    "<img src='https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/images/Pulsars/Geminga_phase_def.png'>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once you have defined your off-pulse interval, use [gtselect](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtselect.txt) to filter out the unwanted events when the pulsar is ON."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "gtselect evclass=128 evtype=3 phasemin=.22 phasemax=.58\n",
    "    ./data/Geminga_gtis.fits\n",
    "    ./data/Geminga_22_58_phasecut.fits\n",
    "    INDEF\n",
    "    INDEF\n",
    "    INDEF\n",
    "    INDEF\n",
    "    INDEF\n",
    "    100\n",
    "    500000\n",
    "    90"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that phasemin and phasemax are not prompted parameters, so you will need to define them on the command line.\n",
    "\n",
    "If your off-pulse interval spans the ends of the phase assignment (as would be the case here if we used the interpulse region) you will need to extract two events segments corresponding to phases 0 to .05, and phases .72 to 1. [gtselect](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtselect.txt) does not handle multiple phase cuts, so you should consider using the FTOOL **ftselect** instead."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 5. Verify the pulsar is missing\n",
    "\n",
    "Now you should bin your data again and check to see if the contribution from the pulsar is really gone. Use [gtbin](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtbin.txt) to bin the remaining data, and view it in *ds9*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "gtbin\n",
    "    CMAP\n",
    "    ./data/Geminga_22_58_phasecut.fits\n",
    "    ./data/Geminga_22_58_phasecut_cmap.fits\n",
    "    ./data/spacecraft.fits\n",
    "    300\n",
    "    300\n",
    "    0.1\n",
    "    CEL\n",
    "    98.476204\n",
    "    17.770661\n",
    "    0\n",
    "    CAR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!ds9 ./data/Geminga_22_58_phasecut_cmap.fits"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src='https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/images/Pulsars/Geminga_OffOn_Compare.png'>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see above, the contribution from Geminga in the new map (left) is mostly gone. The fainter source near the edge of the ROI is now the strongest source in the field. What you will also notice is that you now have many fewer events. That's because the contributions from the other sources has been scaled down by the fraction of data you have removed; in this case the remaining data is 36% (58-22) of the original. This is, in effect, a cut in livetime on these sources."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src='https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/images/Pulsars/Geminga_OffOff_Compare.png'>\n",
    "You can adjust the cuts on phase to try to remove more and more of the pulsar, but this significantly affects the livetime of your other sources, so you should use these cuts with care. Above is a comparison between the previous phase interval of .22-.58 (left) and a tighter cut of .25-.55 (right). As you can see, the pulsar contribution is diminished, but not completely gone."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This data can now be used with a source model in a normal binned or unbinned likelihood analysis to find the flux of the fainter sources (you should still leave a source at the position of Geminga to absorb the residual flux since this pulsar never turns off completely). Once you have fitted your data, multiply the resulting fluxes by 2.778 (or 1/.36) to get the proper flux values for the sources.\n",
    "\n",
    "While Geminga is not really close to other bright gamma-ray sources, regions like Cygnus or Carina contain multiple bright sources that lie very close together. Such regions may benefit significantly from the use of pulsar gating to remove pulsar contributions. However, be aware that every pulsar you remove costs a significant fraction of livetime, so you should use this technique sparingly."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using TEMPO2 for phase gating\n",
    "\n",
    "Since many young gamma-ray pulsars (like Vela) are very noisy, you may need to use TEMPO2, a pulsar analysis package developed by radio astronomers, to assign phases for more complex ephemerides. You can get all the information and download the TEMPO2 software at http://www.atnf.csiro.au/research/pulsar/tempo2/index.php?n=Main.HomePage.\n",
    "\n",
    "A plugin designed to ingest LAT event data and export phase assignments is already included in the TEMPO2 package, so the software should be compatible with your LAT data file. [Instructions](http://fermi.gsfc.nasa.gov/ssc/data/analysis/user/Fermi_plug_doc.pdf) on how to install and run the TEMPO2 package with the LAT plugin have been provided by the developer.\n",
    "\n",
    "For TEMPO2 to run, the TEMPO2 environment variable must be set to the full path to the T2runtime directory. During installation you should specify that the plugins should be installed in this same directory. Otherwise, you will need to copy the plugins directory into this folder."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, let's move the files related to Geminga to some other directory so that they don't get mixed up with Vela:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!mkdir ./data/Geminga\n",
    "!mv ./data/* ./data/Geminga"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we use a recent [Vela parameter file](https://fermi.gsfc.nasa.gov/ssc/data/access/lat/ephems/0835-4510/0835-4510_ApJ_713_154_2010.par) from the [LAT Pulsar Ephemeris](http://fermi.gsfc.nasa.gov/ssc/data/access/lat/ephems/) page."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!wget https://fermi.gsfc.nasa.gov/ssc/data/access/lat/ephems/0835-4510/0835-4510_ApJ_713_154_2010.par\n",
    "!mv *.par ./data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Checking the content of the .par file, we see the ephemeris is valid from MJD 54647.403 to MJD 55014.192.\n",
    "\n",
    "Download the corresponding data from the [LAT Data Server](https://fermi.gsfc.nasa.gov/cgi-bin/ssc/LAT/LATDataQuery.cgi) using these search terms:\n",
    "* Vela (128.287, -45.1901)\n",
    "* 15° ROI\n",
    "* 239557417, 268185600 MET\n",
    "* Energy Range: (100, 500000) MeV"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src='https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/images/Pulsars/Vela_query.png'>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now download the files related to Vela into your current directory.\n",
    "\n",
    "We will assume that the spacecraft file for Vela will be renamed to `VelaSp.fits`. We will also create a list of the photon event files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!mv *.fits ./data\n",
    "!mv ./data/*_SC00.fits ./data/VelaSp.fits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!ls ./data/*_PH*.fits > ./data/Vela.list\n",
    "!cat ./data/Vela.list"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Like Geminga, Vela is very bright so we get several photon files that we need to combine, trim, and apply the zenith cut using [gtselect](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtselect.txt). This step is always followed immediately by [gtmktime](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtmktime.txt):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "gtselect evclass=128 evtype=3\n",
    "    @./data/Vela.list\n",
    "    ./data/Vela_events.fits\n",
    "    INDEF\n",
    "    INDEF\n",
    "    INDEF\n",
    "    INDEF\n",
    "    INDEF\n",
    "    100\n",
    "    500000\n",
    "    90"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "gtmktime\n",
    "    ./data/VelaSp.fits\n",
    "    (DATA_QUAL>0)&&(LAT_CONFIG==1)\n",
    "    no\n",
    "    ./data/Vela_events.fits\n",
    "    ./data/Vela_gtis.fits"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looking at the data, we see the very bright pulsar in the center of the field."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "gtbin\n",
    "    CMAP\n",
    "    ./data/Vela_gtis.fits\n",
    "    ./data/Vela_cmap.fits\n",
    "    ./data/VelaSp.fits\n",
    "    300\n",
    "    300\n",
    "    0.1\n",
    "    CEL\n",
    "    128.836\n",
    "    -45.1763\n",
    "    0\n",
    "    CAR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!ds9 ./data/Vela_cmap.fits"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src='https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/images/Pulsars/Vela_On.png'>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Assigning pulse phases\n",
    "\n",
    "To assign pulse phases, first run TEMPO2 and check if the plugins are detected by:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!tempo2 -h"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You should see \"Fermi\" included in the list of plugins at the end of the help file. To see what options are available with the Fermi plugin, type:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!tempo2 -gr fermi -h"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To assign pulse phases to your FITS file, specify the Fermi plugin, the events (`ft1`) file, the spacecraft (`ft2`) file, and the ephemeris (`-f`).\n",
    "\n",
    "To add a column to the FITS file with the pulse phase, use the `-phase` flag. The `-graph 0` flag suppresses the graphical output from the tool.\n",
    "\n",
    "The Fermi plugin for TEMPO2 was written by Lucas Guillamot, and a [user's guide](http://fermi.gsfc.nasa.gov/ssc/data/analysis/user/Fermi_plug_doc.pdf) is available at the FSSC."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!tempo2 -gr fermi -ft1 ./data/Vela_gtis.fits -ft2 ./data/VelaSp.fits -f 0835-4510_ApJ_713_154_2010.par -phase -graph 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After viewing the lightcurve in *fv*, the off-pulse for Vela seems to be .68-.03. This will require two phase cuts (0-.03 and .68-1). Remember, [gtselect](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtselect.txt) does not handle multiple phase cuts, so we will use **ftselect**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!ftselect ' ./data/Vela_gtis.fits[events]' \\\n",
    "          ./data/Vela_phasecut.fits \\\n",
    "          '(PULSE_PHASE > 0 && PULSE_PHASE < 0.03) || (PULSE_PHASE > 0.68 && PULSE_PHASE < 1)'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you compare the number of events in the pre- and post-filter files, you should find the number of events remaining is about one quarter the original number. Since Vela is the brightest gamma-ray source in the sky, this is to be expected.\n",
    "\n",
    "Now let's see what remains of Vela in our dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "gtbin\n",
    "    CMAP\n",
    "    ./data/Vela_phasecut.fits\n",
    "    ./data/Vela_phasecut_cmap.fits\n",
    "    ./data/spacecraft.fits\n",
    "    300\n",
    "    300\n",
    "    0.1\n",
    "    CEL\n",
    "    128.836\n",
    "    -45.1763\n",
    "    0\n",
    "    CAR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!ds9 ./data/Vela_phasecut_cmap.fits"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src='https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/images/Pulsars/Vela_Off.png'>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By adding another frame and scaling to have the same maximum value, we can compare what the region looks like with and without Vela.\n",
    "\n",
    "<img src='https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/images/Pulsars/Vela_OnOff_Compare.png'>\n",
    "\n",
    "This region contains several high-energy galactic sources that are difficult to detect under the strong galactic diffuse emission that is present. You can solve this by cutting out the low-energy events by using [gtselect](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtselect.txt) and changing the minimum energy to something higher, like 1 GeV."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "gtselect\n",
    "    ./data/Vela_phasecut.fits\n",
    "    ./data/Vela_phasecut_gt1gev.fits\n",
    "    INDEF\n",
    "    INDEF\n",
    "    INDEF\n",
    "    INDEF\n",
    "    INDEF\n",
    "    1000\n",
    "    500000\n",
    "    90"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "gtbin\n",
    "    CMAP\n",
    "    ./data/Vela_phasecut_gt1gev.fits\n",
    "    ./data/Vela_phasecut_gt1gev_cmap.fits\n",
    "    ./data/spacecraft.fits\n",
    "    300\n",
    "    300\n",
    "    0.1\n",
    "    128.836\n",
    "    -45.1763\n",
    "    0\n",
    "    CEL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!ds9 ./data/Vela_phasecut_gt1gev_cmap.fits"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you smooth the image, you will see that the majority of the galactic diffuse emission has been removed.\n",
    "\n",
    "Now, two circular features emerge: VelaX, the pulsar wind nebula powered by the gated pulsar, and Vela Junior, a supernova remnant in the lower left quadrant of the image.\n",
    "\n",
    "<img src='https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/images/Pulsars/Vela_Off_smoothed.png'>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Both these sources lie along the galactic plane and are completely invisible without using gating to remove the strong signal from the pulsar.\n",
    "\n",
    "To measure the significance of these two sources, you will need to perform an [extended source analysis](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/extended/extended.html) to fit the data to a template for each source.\n",
    "\n",
    "Don't forget that you have reduced the content of the data by a significant fraction. The remaining data is only 35% (100-68+3) of the original dataset. You will need to scale the flux results of any likelihood analysis by 100/35=2.857 to recover the flux of the sources from the original dataset."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Proper attribution\n",
    "\n",
    "One final note: The collection and analysis of radio and LAT data that goes into generating each ephemeris is a significant time investment.\n",
    "\n",
    "Please be sure to **acknowledge the source of the ephemeris and any appropriate references** should you use any of these ephemerides for a publication.\n",
    "\n",
    "Please review the LAT Pulsar Ephemerides page for requested acknowledgements, email the timing contact for the pulsar ephemeris you have used, or contact the [Fermi helpdesk](https://fermi.gsfc.nasa.gov/ssc/help/) if you have any questions regarding this policy."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
